Estimativas de Kernel (ou Mapas de Calor)

• Estimador de intensidade de distribuição de pontos:

\[\hat{\lambda}_{\tau}(u) = \dfrac{1}{\tau^2}\sum k(\dfrac{d(u_i , u)}{\tau}), d(u_i , u) \leq \tau\] *Fonte: Referência científica: Druck, S.; Carvalho, M.S.; Câmara, G.; Monteiro, A.V.M. (eds) “Análise Espacial de Dados Geográficos”. Brasília, EMBRAPA, 2004

• Localização da ocorrência de casos de dengue em Belo Horizonte/MG

Mapa Temático

• Taxas de câncer de pulmão na população branca masculina nos Estados Unidos, por condados no ano de 1998

Matriz de Vizinhança

Moran Global, Moran Local e Lisa Map

Moran Global

Moran Local

• Desigualdade no nível distrital na cobertura de saúde reprodutiva, materna, neonatal e infantil na Índia

Fonte: PANDA, Basant Kumar; KUMAR, Gulshan; AWASTHI, Ashish. District level inequality in reproductive, maternal, neonatal and child health coverage in India. BMC public health, v. 20, n. 1, p. 1-10, 2020.

Modelos de Regressão Espacial

• Hipótese de independência das observações em geral é Falsa \(\rightarrow\) Dependência Espacial

• Efeitos Espaciais \(\rightarrow\) Se existir forte tendência ou correlação espacial, os resultados serão influenciados, apresentando associação estatística onde não existe (e vice-versa);

• Como verificar ? \(\rightarrow\) Medir a autocorrelação espacial dos resíduos da regressão (Índice de Moran dos resíduos).

• Autocorrelação espacial constatada ! E agora ? \(\rightarrow\) Utilizar Modelos de regressão que incorporam efeitos espaciais:

  1. Modelos Globais: utilizam um único parâmetro para capturar a estrutura de correlação espacial \(\rightarrow\) ex: Spatial Error Models (CAR)

\[y_i = \beta_0 + \sum^{p}_{k} \beta_k x_{ik} + \varepsilon_i \text{ , sendo } \varepsilon_i = \lambda W + \xi\]

• Os efeitos da autocorrelação espacial são associados ao termo de erro. Sendo \(W\) a matriz de vizinhaça, \(\varepsilon\) a componente do erro com efeitos espaciais, \(λ\) é o coeficiente autoregressivo e \(\xi\) é a componente do erro com variância constante e não correlacionada.

• A hipótese nula para a não existência de autocorrelação é que \(\lambda = 0\), ou seja, o termo de erro não é espacialmente correlacionado.

  2. Modelos Locais: parâmetros variam continuamente no espaço ex: Geographically Weighted Regression (GWR)

\[y_i = \beta_{0}(u_i, v_i) + \sum^{p}_{k} \beta_k (u_i, v_i) x_{ik} + \varepsilon_i \text{ , sendo } (u_i, v_i) \text{ as coordenadas geográficas}\] Fonte: ARDILLY, Pascal et al. Manuel d’analyse spatiale. 2018.

SiG QGIS

QGIS: Um Sistema de Informação Geográfica livre e aberto

GEODA

GEODA: AN INTRODUCTION TO SPATIAL DATA ANALYSIS

R

Fonte: Dicas para integração e instalação do R 4.2 no Ubuntu 22.04 LTS e os pacotes espaciais

Python

Fonte: Github Geopandas

Bibliotecas que serão utilizadas:

library(ggplot2)
library(dplyr)
library(viridis)
library(geobr)
library(sf)
library(maptools)
library(leaflet)
library(ggspatial)

• Para acessar os dados dos limites territoriais de todos os estados brasileiros é necessário utilizar a função read_state.

brasil.ufs <- read_state(code_state = "all", year = 2019, showProgress = FALSE)

Primeiro, vamos fazer um gráfico apenas com as geometrias.

ggplot(brasil.ufs) + geom_sf()

• Para a construção de um mapa onde cada estado recebe uma cor de acordo com a sua região geogrática, procedemos da seguinte forma:

ggplot(brasil.ufs) + geom_sf(aes(fill = name_region)) + labs(fill = "Região")

• Agora utilizaremos dados relativos a porcentagem de municípios com rede de esgoto de acordo com a unidade da federação (fonte dos dados ) para fazer um gráfico. Vamos associar esses dados a tabela de acordo com a variável State e padronizaremos a porcentagem para variar de 0 a 100.

# Entrando com os dados observados na wikipedia
acesso_san <- data.frame(code_state = c(12, 27, 16, 13, 29, 23, 53, 32, 52, 21, 51,
    50, 31, 15, 25, 41, 26, 22, 33, 24, 43, 11, 14, 42, 35, 28, 17), com_rede = c(0.273,
    0.412, 0.313, 0.177, 0.513, 0.696, 1, 0.974, 0.28, 0.065, 0.191, 0.449, 0.916,
    0.063, 0.731, 0.421, 0.881, 0.045, 0.924, 0.353, 0.405, 0.096, 0.4, 0.352, 0.998,
    0.347, 0.129))

# Construindo o mapa com ggplot
brasil.ufs |>
    left_join(acesso_san, by = "code_state") |>
    mutate(com_rede = 100 * com_rede) |>
    ggplot(aes(fill = com_rede), color = "black") + geom_sf() + scale_fill_viridis(name = "Municípios com rede de esgoto (%)",
    direction = -1) + xlab("Longitude") + ylab("Latitude") + annotation_scale(location = "bl") +
    annotation_north_arrow(location = "br") + theme_minimal()

• Uma forma alteranativa de apresentar esses mesmos dados se dá pela apresentação de círculos com raios proporcionais a porcentgem de municípios com rede de esgoto no centroide de cada geometria (nesse caso, UF).

# Juntando o banco com os dados + malha
coord_pontos <- brasil.ufs |>
    left_join(acesso_san, by = "code_state") |>
    mutate(com_rede = 100 * com_rede) |>
    st_centroid()

# Construindo o mapa com ggplot
ggplot(brasil.ufs) + geom_sf() + geom_sf(data = coord_pontos, aes(size = com_rede),
    col = "blue", alpha = 0.65, show.legend = "point") + scale_size_continuous(name = "Municípios com rede de esgoto (%)") +
    xlab("Longitude") + ylab("Latitude") + annotation_scale(location = "bl") + annotation_north_arrow(location = "br") +
    theme_minimal()

• Uma alternativa interativa para trabalhar com mapas é com a utilização do pacote leaflet

data.frame(st_coordinates(coord_pontos), com_rede = coord_pontos$com_rede, UF = coord_pontos$name_state) |>
    leaflet() |>
    addTiles() |>
    addCircleMarkers(~X, ~Y, label = ~as.character(paste0(UF, ": ", com_rede, "%")),
        labelOptions = labelOptions(textsize = "13px"), radius = ~com_rede/10, fillOpacity = 0.5)

Créditos ao Thiago Mendonça

Biliotecas do R que serão utilizadas

library(readr)
library(tidyverse)
library(sf)
library(maptools)
library(spatstat)
library(tmap)

Lendo a tabela da população por setor censitário e os shapes files do contorno e por setor censitário de Dourados/MS

pop2010 <- read_csv("dados/dengue_dourados/pop2010.csv")
setor.sf <- read_sf("dados/dengue_dourados/Setor_UTM_SIRGAS.shp", crs = 31981)
contorno.sf <- read_sf("dados/dengue_dourados/contorno.shp", crs = 31981)

OBS: Um aspecto muito importante dos dados espaciais é o sistema de referência de coordenadas (CRS) que é usado. Por exemplo, uma localização de (140, 12) não é significativa se você sabe onde está a origem e se a coordenada x está a 140 metros, quilômetros ou talvez graus de distância dela (na direção x). (link)

Fazendo um join com as tabelas com os setores censitários + população

setor.sf <- setor.sf %>%
    mutate(idsetor = as.numeric(CD_GEOCODI)) %>%
    inner_join(pop2010, by = "idsetor")

Lendo e plotando os casos de dengue georreferenciados em Dourados/MS

casos <- read_csv("dados/dengue_dourados/dengue_dourados.csv")
casos.sf <- st_as_sf(casos, coords = c("X", "Y"), crs = 31981)

ggplot(setor.sf) + geom_sf(fill = "white", color = "black") + geom_sf(data = casos.sf,
    color = "red", size = 1) + theme_void()

Lendo e plotando os os pontos de coleta de lixo georreferenciados em Dourados/MS

lixo <- read_csv2("dados/dengue_dourados/lixo_dourados.csv")
lixo.sf <- st_as_sf(lixo, coords = c("X", "Y"), crs = 31981)

ggplot(setor.sf) + geom_sf(fill = "white", color = "black") + geom_sf(data = lixo.sf,
    color = "blue", size = 1) + theme_void()

• Como podemos observar, existem alguns pontos de coleta de lixo fora do contorno de Dourados/MS

Uma forma de ficarmos só com os pontos dentro do polígono é utilizando o comando st_intersection.

lixo2.sf <- st_intersection(contorno.sf, lixo.sf)

ggplot(setor.sf) + geom_sf(fill = "white", color = "black") + geom_sf(data = lixo2.sf,
    color = "blue", size = 1) + theme_void()

Utilizando as informações dos casos (pontos) + do lixo (ponto) + população de cada setor censitário (mapa temático)

ggplot(setor.sf) + geom_sf(aes(fill = pop)) + geom_sf(data = casos.sf, color = "red",
    size = 0.7, aes(colour = "Caso"), show.legend = "point") + geom_sf(data = lixo2.sf,
    color = "salmon", size = 1, aes(colour = "Lixo"), show.legend = "point") + scale_fill_distiller(palette = "PuBu",
    direction = 1) + scale_colour_manual(values = c(Caso = "red", Lixo = "slamon")) +
    theme_minimal()

#### Construindo os buffers

• Agora serão construidos buffers de 500m de distância ao redor de cada ponto de coleta de lixo. Isso é interessante para verificar se os casos de dengue ocorrem em um raio de até 500m de distância dos pontos de coleta de lixo. Ou seja, a pergunta é, será que a distância dos pontos de coeta de lixo influenciam a ocorrência do caso de dengue ?

Buffers: São polígonos que contornam um objeto a uma determinada distância. Sua principal função é materializar os conceitos de “perto” e “longe”.

lixo_buffer <- st_buffer(lixo2.sf, 500)

ggplot(setor.sf) + geom_sf(aes(fill = pop)) + geom_sf(data = lixo_buffer, color = "gray",
    fill = "transparent", size = 0.4) + geom_sf(data = casos.sf, color = "red", size = 0.7) +
    scale_fill_distiller(palette = "PuBu", direction = 1) + scale_colour_manual(values = c(Caso = "red",
    Lixo = "slamon")) + theme_minimal()

Represntando os casos e o lixo de forma interativa

tm_shape(setor.sf) + tm_borders("black") + tm_shape(casos.sf) + tm_dots("red") +
    tm_shape(lixo.sf) + tm_dots("green") + tm_shape(lixo_buffer) + tm_borders("blue") +
    tmap_mode("view")

Convertendo o dado de pontos (padrão pontual) para dados de área

## conta casos por setor 
casos.sf$contador <- 1 
casos <- setor.sf |>  
  st_join(casos.sf) |> 
  filter(CLASSI_FIN == 1) %>%  ## seleciona somente os casos confirmados
  group_by(ID) |> 
  summarise(casos = sum(contador))

st_geometry(casos) <- NULL  ## remove atributos de geometria

## numero de depositos de Lixo por setor 
lixo.sf$contador <- 1 

lixo <- setor.sf |>  
  st_join(lixo.sf) |> 
  group_by(ID) |> 
  summarise(lixo = sum(contador))

st_geometry(lixo) <- NULL ## remove atributos de geometria

# Inserindo as contagens dos casos e de pontos de coleta de lixo no atributo com a geometria.
setor.sf <- setor.sf |> 
  left_join(lixo,by = 'ID') |> 
  left_join(casos,by = 'ID') 

Plotando o mapa temático dos casos por setor censitário

plot(setor.sf["casos"])

Plotando o mapa temático dos pontos de coleta de lixo por setor censitário

plot(setor.sf["lixo"])

Calculando a taxa de incidência e plotando o mapa temático dos pontos de coleta de lixo por setor censitário

setor.sf$tx <- (setor.sf$casos/setor.sf$pop) * 1000
setor.sf$tx[is.na(setor.sf$tx)] <- 0  # Transformando os missings em zero

summary(setor.sf$tx)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   0.000   1.736   4.065   4.311  56.399 

Plotando a distribuição da incidência em Dourados/MS

library(wesanderson)
pal <- wes_palette("Moonrise3", 20, type = "continuous")

ggplot(setor.sf) + geom_sf(aes(fill = tx), color = "black") + scale_fill_gradientn(colours = pal) +
    ggtitle("Taxa de incidência de Dengue") + theme_void()

Kernel por atributos

• Vamos plotar o kernel por atributos referente a taxa de incidência de dengue em Dourados/MS.

• Primeiramente é necessário dissolver os polígonos em formato sf para obter o contorno. Nesse caso queremos preservar o atributo AREA

dourados.contorno <- st_union(setor.sf)
plot(dourados.contorno)

dourados.w <- as.owin(st_geometry(dourados.contorno))

Extraindo os centróides dos polígonos em Dourados/MS

centroides <- st_centroid(st_geometry(setor.sf))

# Transformando em os centróides em formato sp
centroides.sp <- as.data.frame(as_Spatial(centroides))
names(centroides.sp) <- c("X", "Y")

plot(centroides.sp)

Colocando os pontos no formato sp

centroides.ppp <- ppp(centroides.sp$X, centroides.sp$Y, dourados.w)

plot(centroides.ppp, pch = 19, cex = 0.5)

Fazendo o kernel por atributo da taxa de detecção

kernel.tx <- density(centroides.ppp, 500, weights = setor.sf$tx, scalekernel = TRUE)
plot(kernel.tx)

Construindo a matriz de vizinhança para verificar a autocorrelação espacial

library(spdep)
viz <- poly2nb(setor.sf)
viz

Neighbour list object:
Number of regions: 284 
Number of nonzero links: 1726 
Percentage nonzero weights: 2.139952 
Average number of links: 6.077465 

• Iremos precisar da coordenadas dos centróides

setor.sp <- as(setor.sf, "Spatial")  # convertendo em formato sp
coord <- coordinates(setor.sp)  # coordenadas dos centroidas dos poligonos de dourados
class(setor.sp)

[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"

Verificando a malha de conectividade da vizinhança de Dourados/MS

viz.sf <- as(nb2lines(viz, coords = coord), "sf")
viz.sf <- st_set_crs(viz.sf, st_crs(setor.sf))

# Plota o grafo de conectividade por contiguidade
mapa.viz <- ggplot(setor.sf) + geom_sf(fill = "salmon", color = "white") + geom_sf(data = viz.sf) +
    theme_minimal() + ggtitle("Vizinhança por \n conectividade") + ylab("Latitude") +
    xlab("Longitude")
mapa.viz

Obtendo a correlação da taxa de incidência de dengue Dourados/MS

pesos.viz <- nb2listw(viz)
moran.test(setor.sf$tx, pesos.viz)

    Moran I test under randomisation

data:  setor.sf$tx  
weights: pesos.viz    

Moran I statistic standard deviate = 15.724, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.524720303      -0.003533569       0.001128660 

Plotando o correlograma

correl <- sp.correlogram(viz, setor.sf$tx, order = 8, method = "I")
correl

Spatial correlogram for setor.sf$tx 
method: Moran's I
           estimate expectation    variance standard deviate Pr(I) two sided
1 (284)  0.52472030 -0.00353357  0.00112866          15.7239       < 2.2e-16
2 (284)  0.23776816 -0.00353357  0.00048540          10.9525       < 2.2e-16
3 (284)  0.09536593 -0.00353357  0.00028358           5.8729       4.282e-09
4 (284) -0.02063378 -0.00353357  0.00019736          -1.2172         0.22351
5 (284) -0.07589158 -0.00353357  0.00016732          -5.5939       2.221e-08
6 (284) -0.06448448 -0.00353357  0.00016165          -4.7940       1.635e-06
7 (284) -0.02708340 -0.00353357  0.00016725          -1.8210         0.06861
8 (284) -0.02744543 -0.00353357  0.00018328          -1.7662         0.07735
           
1 (284) ***
2 (284) ***
3 (284) ***
4 (284)    
5 (284) ***
6 (284) ***
7 (284) .  
8 (284) .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(correl)

Mapeando os polígonos que tiveram os p-valores mais significativos no Moran Local.

setor.sf$pval <- localmoran(setor.sf$tx, pesos.viz)[, 5]

tm_shape(setor.sf) + tm_polygons(col = "pval", title = "p-valores", breaks = c(0,
    0.01, 0.05, 0.1, 1), border.col = "white", palette = "-Oranges") + tm_scale_bar(width = 0.15)

Moran Local (Lisa Map) da taxa de incidência Dourados/MS

resI <- localmoran.sad(lm(setor.sf$tx ~ 1), 1:length(viz), viz, style = "W")
summary(resI)[1:10, ]

    Local Morans I Stand. dev. (N)   Pr. (N) Saddlepoint Pr. (Sad)  Expectation
1 1   -0.009885292     -0.01575255 1.0125682 -0.03600714 0.9712767 -0.003533569
2 2    0.226981101      0.46511580 0.6418485  0.70056077 0.4835772 -0.003533569
3 3   -0.010910944     -0.01829621 1.0145974 -0.04041673 0.9677609 -0.003533569
4 4    0.484675519      1.31014141 0.1901480  1.43930439 0.1500643 -0.003533569
5 5    0.018648578      0.05952726 0.9525322  0.09384037 0.9252360 -0.003533569
     Variance    Skewness Kurtosis   Minimum  Maximum         omega       sad.r
1 1 0.1625854 -0.05147857 8.753481 -57.75455 56.75455 -0.0001369880 -0.01569409
2 2 0.2456263 -0.04188294 8.752890 -70.87400 69.87400  0.0028119325  0.44472643
3 3 0.1625854 -0.05147857 8.753481 -57.75455 56.75455 -0.0001590844 -0.01822753
4 4 0.1388594 -0.05570264 8.753780 -53.41199 52.41199  0.0066304485  1.08297547
5 5 0.1388594 -0.05570264 8.753780 -53.41199 52.41199  0.0005595065  0.05929966
          sad.u
1 1 -0.01569909
2 2  0.49831660
3 3 -0.01823490
4 4  1.59298211
5 5  0.05942125
 [ reached 'max' / getOption("max.print") -- omitted 5 rows ]

setor.sf$MoranLocal <- summary(resI)[, 1]

library(scales)

ggplot(setor.sf) + geom_sf(aes(fill = MoranLocal), color = "black") + scale_fill_gradientn(colours = c("blue",
    "white", "red"), values = rescale(c(min(setor.sf$MoranLocal), 0, max(setor.sf$MoranLocal))),
    guide = "colorbar") + ggtitle("Moran local") + theme_void()

Ajustando o modelo de regressão linear simples.

# Transformando os missings em zero
setor.sf$lixo[is.na(setor.sf$lixo)] <- 0

dourados.lm <- lm(tx ~ lixo, data = setor.sf)
summary(dourados.lm)

Call:
lm(formula = tx ~ lixo, data = setor.sf)

Residuals:
   Min     1Q Median     3Q    Max 
-4.510 -4.115 -2.286  0.237 51.889 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   4.5103     0.4889   9.225   <2e-16 ***
lixo         -0.3955     0.1945  -2.033    0.043 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.365 on 282 degrees of freedom
Multiple R-squared:  0.01445,   Adjusted R-squared:  0.01095 
F-statistic: 4.134 on 1 and 282 DF,  p-value: 0.04298

Checando os residuos para verificar a presença de autocorrelação.

dourados.lm$lmresid <- residuals(dourados.lm)
moran.test(dourados.lm$lmresid, pesos.viz)

    Moran I test under randomisation

data:  dourados.lm$lmresid  
weights: pesos.viz    

Moran I statistic standard deviate = 15.578, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.520064269      -0.003533569       0.001129669 

Ajustando o modelo CAR (Spatial Error Model)

library(spatialreg)
dourados.car <- errorsarlm(tx ~ lixo, data = setor.sf, listw = pesos.viz)
summary(dourados.car)

Call:errorsarlm(formula = tx ~ lixo, data = setor.sf, listw = pesos.viz)

Residuals:
      Min        1Q    Median        3Q       Max 
-12.53213  -2.34032  -1.09686   0.82021  31.62992 

Type: error 
Coefficients: (asymptotic standard errors) 
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  4.07199    1.54998  2.6271 0.008611
lixo        -0.25129    0.14814 -1.6963 0.089826

Lambda: 0.8063, LR test value: 167.98, p-value: < 2.22e-16
Asymptotic standard error: 0.041575
    z-value: 19.394, p-value: < 2.22e-16
Wald statistic: 376.13, p-value: < 2.22e-16

Log likelihood: -885.0643 for error model
ML residual variance (sigma squared): 25.435, (sigma: 5.0433)
Number of observations: 284 
Number of parameters estimated: 4 
AIC: 1778.1, (AIC for lm: 1944.1)

Checando os residuos para verificar a presença de autocorrelação

dourados.car$carresid <- residuals(dourados.car)
moran.test(dourados.car$carresid, pesos.viz)

    Moran I test under randomisation

data:  dourados.car$carresid  
weights: pesos.viz    

Moran I statistic standard deviate = 0.78357, p-value = 0.2166
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.023005067      -0.003533569       0.001147099 

Ajustando o modelo GWR (Geographically Weighted Regression)

# Biblioteca para ajustar o modelos GWR
library(spgwr)
# Estimando a largura de banda “ideal” para o kernel
GWRbanda <- gwr.sel(tx ~ lixo, data = setor.sf, coords = cbind(centroides.sp$X, centroides.sp$Y),
    adapt = T)

Adaptive q: 0.381966 CV score: 14699.46 
Adaptive q: 0.618034 CV score: 15071.4 
Adaptive q: 0.236068 CV score: 14204.85 
Adaptive q: 0.145898 CV score: 13336.85 
Adaptive q: 0.09016994 CV score: 12151.44 
Adaptive q: 0.05572809 CV score: 10995.75 
Adaptive q: 0.03444185 CV score: 10093.37 
Adaptive q: 0.02128624 CV score: 9758.837 
Adaptive q: 0.01315562 CV score: 9295.706 
Adaptive q: 0.008130619 CV score: 8996.895 
Adaptive q: 0.005024999 CV score: 8988.675 
Adaptive q: 0.00638842 CV score: 9000.179 
Adaptive q: 0.00310562 CV score: 9842.835 
Adaptive q: 0.004291861 CV score: 9067.819 
Adaptive q: 0.005545779 CV score: 8976.568 
Adaptive q: 0.005867639 CV score: 8980.638 
Adaptive q: 0.005586469 CV score: 8976.67 
Adaptive q: 0.005505089 CV score: 8976.598 
Adaptive q: 0.005545779 CV score: 8976.568 

# Ajustando o modelo GWR
dourados.gwr = gwr(tx ~ lixo, data = setor.sf, coords = cbind(centroides.sp$X, centroides.sp$Y),
    adapt = GWRbanda, hatmatrix = TRUE, se.fit = TRUE)

dourados.gwr

Call:
gwr(formula = tx ~ lixo, data = setor.sf, coords = cbind(centroides.sp$X, 
    centroides.sp$Y), adapt = GWRbanda, hatmatrix = TRUE, se.fit = TRUE)
Kernel function: gwr.Gauss 
Adaptive quantile: 0.005545779 (about 1 of 284 data points)
Summary of GWR coefficient estimates at data points:
                   Min.    1st Qu.     Median    3rd Qu.       Max.  Global
X.Intercept.  -0.324952   1.412188   3.099163   5.397652  37.581278  4.5103
lixo         -18.321709  -1.222217  -0.400929  -0.061492   1.555744 -0.3955
Number of data points: 284 
Effective number of parameters (residual: 2traceS - traceS'S): 134.4484 
Effective degrees of freedom (residual: 2traceS - traceS'S): 149.5516 
Sigma (residual: 2traceS - traceS'S): 5.447405 
Effective number of parameters (model: traceS): 100.4864 
Effective degrees of freedom (model: traceS): 183.5136 
Sigma (model: traceS): 4.917575 
Sigma (ML): 3.952993 
AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 1904.233 
AIC (GWR p. 96, eq. 4.22): 1687.144 
Residual sum of squares: 4437.827 
Quasi-global R2: 0.7140881 

# Colocando a saída do modelo dentro de um objeto dataframe.
results <- as.data.frame(dourados.gwr$SDF)
head(results)

     sum.w X.Intercept.       lixo X.Intercept._se  lixo_se       gwr.e
1 4.630090     6.489949 -1.8954895        2.411594 1.718037 -0.02159391
2 3.934954     7.883872 -1.2811203        2.504116 1.630176  1.98454873
3 3.062602     7.267798 -1.5625806        2.780695 1.625888 -1.85906371
4 4.139484     9.791258 -1.1687314        2.376803 1.277302  6.80957181
5 5.541437     4.586352 -0.6839813        2.090863 1.185501 -2.63010618
      pred  pred.se   localR2 X.Intercept._se_EDF lixo_se_EDF pred.se.1
1 2.698970 2.477284 0.6403656            2.671425    1.903141  2.744191
2 7.883872 2.504116 0.5043081            2.773914    1.805815  2.773914
3 5.705218 2.164469 0.5288378            3.080292    1.801064  2.397673
4 8.622527 1.794393 0.4101644            2.632885    1.414921  1.987725
5 3.902371 1.524868 0.6665163            2.316137    1.313229  1.689160
   coord.x coord.y
1 731406.5 7541547
2 730845.4 7541481
3 730957.7 7541239
4 730487.2 7541437
5 729874.1 7541446
 [ reached 'max' / getOption("max.print") -- omitted 1 rows ]

Verificando a distribuição dos coeficientes de regressão para a variável lixo

hist(results$lixo)
abline(v = median(results$lixo), col = "red")

Verificando a distribuição dos localR2

hist(results$localR2)
abline(v = median(results$localR2), col = "blue")

Incorporando alguns parâmetros de saída do modelo na tabela setor.sf

setor.sf$coef.lixo <- results$lixo
setor.sf$localR2 <- results$localR2
setor.sf$pred.gwr <- results$pred

Mapa dos coeficientes de regressão para a variável lixo

map.lixo <- ggplot(setor.sf) + geom_sf(aes(fill = coef.lixo), color = "black") +
    scale_fill_gradientn(colours = pal) + ggtitle("Distribuição dos coef. var. lixo") +
    theme_void()
map.lixo

Checando os residuos para verificar a presença de autocorrelação para o modelo GWR.

# Calculando os resíduos para o modelo GWR
results$residuos <- setor.sf$tx - results$pred

moran.test(results$residuos, pesos.viz)

    Moran I test under randomisation

data:  results$residuos  
weights: pesos.viz    

Moran I statistic standard deviate = 1.4806, p-value = 0.06935
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.046816835      -0.003533569       0.001156432 

Mapeando os coeficientes de regressão para a variável lixo por significancia através do teste de wald

# Calculando a estatística de wald
setor.sf$wald.teste <- abs(results$lixo/results$lixo_se)
# Dicotomizando a estatística de wald
setor.sf$wald.teste <- ifelse(setor.sf$wald.teste < 2, 0, 1)


map.wald <- ggplot(setor.sf) + geom_sf(aes(fill = factor(wald.teste)), color = "black") +
    scale_fill_manual(values = c("white", "purple"), labels = c("< 2", ">= 2"), name = "Wald") +
    ggtitle("Coef. lixo significativos") + theme_void()


library(gridExtra)
grid.arrange(map.lixo, map.wald, ncol = 2)

Mapa dos coeficientes de determinação regionalizados (\(R^2\) local).

ggplot(setor.sf) + geom_sf(aes(fill = localR2), color = "black") + scale_fill_gradientn(colours = pal) +
    ggtitle("R² local") + theme_void()

Verificando a distribuição dos preditos.

histdens <- function(x, titulo = "") {
    densi <- density(x)
    xli <- range(densi$x)
    yli <- range(densi$y)
    hist(x, col = "red", probability = T, xlim = xli, ylim = yli, main = titulo)
    lines(densi, lwd = 2)
    abline(v = median(x), lwd = 4, col = 4, lty = 2)
}

par(mfrow = c(2, 2))

hist.tx <- histdens(setor.sf$tx, "Tx Bruta")
hist.lm <- histdens(dourados.lm$fitted.values, "Pred LM")
hist.car <- histdens(dourados.car$fitted.values, "Pred CAR")
hist.gwr <- histdens(results$pred, "Pred GWR")

Mapeando os valores observados e preditos dos modelos ajustados

library(colorspace)  # 

setor.sf$brks <- cut(setor.sf$tx, include.lowest = TRUE, right = TRUE, breaks = c(-4,
    0, 2, 4, 10, 57), labels = c("0", "0 - 2", "2 - 4", "4 - 10", "> 10"))

tx.bruta <- ggplot(setor.sf) + geom_sf(aes(fill = brks), color = "black") + ggtitle("Taxa Bruta") +
    scale_fill_discrete_sequential(palette = "Heat2", c1 = 80, c2 = 30, l1 = 30,
        l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75", drop = FALSE, name = "Taxa") +
    theme_void()



setor.sf$brks.lm <- cut(dourados.lm$fitted.values, lowest = TRUE, right = TRUE, breaks = c(-4,
    0, 2, 4, 10, 57), labels = c("0", "0 - 2", "2 - 4", "4 - 10", "> 10"))


pred.lm <- ggplot(setor.sf) + geom_sf(aes(fill = brks.lm), color = "black") + ggtitle("Taxa Predita - LM") +
    scale_fill_discrete_sequential(palette = "Heat2", c1 = 80, c2 = 30, l1 = 30,
        l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75", drop = FALSE, name = "Taxa") +
    theme_void()

setor.sf$brks.car <- cut(dourados.car$fitted.values, lowest = TRUE, right = TRUE,
    breaks = c(-4, 0, 2, 4, 10, 57), labels = c("0", "0 - 2", "2 - 4", "4 - 10",
        "> 10"))


pred.car <- ggplot(setor.sf) + geom_sf(aes(fill = brks.car), color = "black") + ggtitle("Taxa Predita - CAR") +
    scale_fill_discrete_sequential(palette = "Heat2", c1 = 80, c2 = 30, l1 = 30,
        l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75", drop = FALSE, name = "Taxa") +
    theme_void()

setor.sf$brks.gwr <- cut(results$pred, lowest = TRUE, right = TRUE, breaks = c(-4,
    0, 2, 4, 10, 57), labels = c("0", "0 - 2", "2 - 4", "4 - 10", "> 10"))


pred.gwr <- ggplot(setor.sf) + geom_sf(aes(fill = brks.car), color = "black") + ggtitle("Taxa Predita - GWR") +
    scale_fill_discrete_sequential(palette = "Heat2", c1 = 80, c2 = 30, l1 = 30,
        l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75", drop = FALSE, name = "Taxa") +
    theme_void()

library(gridExtra)
grid.arrange(tx.bruta, pred.lm, pred.car, pred.gwr, ncol = 2)

Comparando a distribuição dos resíduos dos modelos ajustados

library(vioplot)
vioplot(dourados.lm$residuals, dourados.car$residuals, results$residuos, names = c("LM",
    "CAR", "GWR"), col = "orange")
title("Gráficos violinos da distribuição dos resíduos")

Mapeando os resíduos dos modelos ajustados

library(colorspace)  # 

setor.sf$brks.res.lm <- cut(dourados.lm$residuals, include.lowest = TRUE, right = TRUE,
    breaks = c(-14, -5, -1, 1, 5, 52), labels = c("< -5", "-5 a -1", "0", "1 a 5",
        "> 5"))

res.lm <- ggplot(setor.sf) + geom_sf(aes(fill = brks.res.lm), color = "black") +
    ggtitle("Resíduos - LM") + scale_fill_discrete_sequential(palette = "Purple-Yellow",
    c1 = 80, c2 = 30, l1 = 30, l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75",
    drop = FALSE, name = "Taxa") + theme_void()

setor.sf$brks.res.car <- cut(dourados.car$residuals, include.lowest = TRUE, right = TRUE,
    breaks = c(-14, -5, -1, 1, 5, 52), labels = c("< -5", "-5 a -1", "0", "1 a 5",
        "> 5"))

res.car <- ggplot(setor.sf) + geom_sf(aes(fill = brks.res.car), color = "black") +
    ggtitle("Resíduos - CAR") + scale_fill_discrete_sequential(palette = "Purple-Yellow",
    c1 = 80, c2 = 30, l1 = 30, l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75",
    drop = FALSE, name = "Taxa") + theme_void()


setor.sf$brks.res.gwr <- cut(results$residuos, include.lowest = TRUE, right = TRUE,
    breaks = c(-14, -5, -1, 1, 5, 52), labels = c("< -5", "-5 a -1", "0", "1 a 5",
        "> 5"))

res.gwr <- ggplot(setor.sf) + geom_sf(aes(fill = brks.res.gwr), color = "black") +
    ggtitle("Resíduos - GWR") + scale_fill_discrete_sequential(palette = "Purple-Yellow",
    c1 = 80, c2 = 30, l1 = 30, l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75",
    drop = FALSE, name = "Taxa") + theme_void()


library(gridExtra)
grid.arrange(res.lm, res.car, res.gwr, ncol = 2)

Biliotecas do R que serão utilizadas

library(sf)
library(sp)
library(dplyr)
library(gstat)
library(lattice)
library(automap)
library(raster)
library(leaflet)

Importando a tabela com a chuva acumulada média de 7 dias das últimas 24hs e 96hs das estações pluviométricas da cidde do Rio de Janeiro

Descrição das Estações (Alerta Rio)

Download Dados Pluviométricos (Alerta Rio)

pluvio <- read.csv("dados/chuva_rio/pluviosidade.csv", sep = ";")

Análise Gráfica Descritiva

quantil <- quantile(pluvio$acumulado_24h, seq(0, 1, 0.2))
quantil

   0%   20%   40%   60%   80%  100% 
 0.00  1.92  4.24  7.36 11.84 36.00 

Transformar os dados em um objeto espacial do R

• x - Longitude
• y - Latitude

sp::coordinates(pluvio) = ~x + y

Análise Gráfica Descritiva com os dados espaciais

# Bubble plot
bubble(pluvio, "acumulado_24h", key.entries = quantil, pch = 19, col = "blue")

# Point plot
spplot(pluvio["acumulado_24h"], scales = list(draw = T), key.space = "right", colorkey = T)

Modelando o variograma experimental (ou empírico)

• width - Distância média entre amostras ou distância dos lags
• cutoff - Máxima distância

variogram.emp = variogram(acumulado_24h ~ x + y, pluvio, width = 1000, cutoff = 20000)
variogram.emp

   np      dist     gamma dir.hor dir.ver   id
1   2  1385.875 20.381155       0       0 var1
2   6  2543.894  8.139327       0       0 var1
3   8  3648.158  5.261490       0       0 var1
4  15  4583.562 13.707344       0       0 var1
5  22  5491.324 26.228076       0       0 var1
6  15  6531.432 52.574023       0       0 var1
7  12  7512.398 42.179960       0       0 var1
8  20  8516.447 65.283855       0       0 var1
9  15  9432.920 36.022491       0       0 var1
10 19 10450.957 36.569808       0       0 var1
11 13 11492.778 76.150012       0       0 var1
12 19 12568.517 39.354425       0       0 var1
 [ reached 'max' / getOption("max.print") -- omitted 7 rows ]

# Variogram plot
plot(variogram.emp, main = "Empirical variogram", pch = 19, col = "darkblue")

Ajustando o semivariograma teórico

• Sill - Semivariância estrutural ou contribuição (ponto máximo que chega ao plato no eixo de y)
• Range - Alcance (ponto máximo em x)
• Nugget - Efeito pepita, quando (y, 0)

Fonte do gráfico

## 1) Modelo Linear
lin.fit = fit.variogram(variogram.emp, model = vgm(psill = 240, model = "Lin", range = 5000,
    nugget = 10))

## 2) Modelo exponencial
exp.fit = fit.variogram(variogram.emp, model = vgm(psill = 240, model = "Exp", range = 5000,
    nugget = 10))

## 3) Modelo gaussiano
gau.fit = fit.variogram(variogram.emp, model = vgm(psill = 240, model = "Gau", range = 5000,
    nugget = 10))

## 3) Modelo wave
wav.fit = fit.variogram(variogram.emp, model = vgm(psill = 240, model = "Wav", range = 5000,
    nugget = 10))

plot.lin <- plot(variogram.emp, lin.fit, main = "Modelo Linear")
plot.exp <- plot(variogram.emp, exp.fit, main = "Modelo Exponencial")
plot.gau <- plot(variogram.emp, gau.fit, main = "Modelo Gaussiano")
plot.wav <- plot(variogram.emp, wav.fit, main = "Modelo Wave")

Validação cruzada

• RMSE (root mean squared error): é a medida que calcula “a raiz quadrática média” dos erros entre valores observados (reais) e predições (hipóteses).

## 1) Modelo Linear
cv.lin <- krige.cv(acumulado_24h ~ x + y, locations = pluvio, model = lin.fit)
summary(cv.lin)

Object of class SpatialPointsDataFrame
Coordinates:
        min     max
x  634918.7  687966
y 7450222.8 7475508
Is projected: NA 
proj4string : [NA]
Number of points: 30
Data attributes:
   var1.pred          var1.var         observed        residual        
 Min.   : 0.6757   Min.   : 9.001   Min.   : 0.00   Min.   :-10.41265  
 1st Qu.: 4.3063   1st Qu.:11.956   1st Qu.: 2.20   1st Qu.: -2.84581  
 Median : 6.7569   Median :19.120   Median : 5.00   Median : -0.44588  
 Mean   : 7.2750   Mean   :28.832   Mean   : 7.34   Mean   :  0.06504  
 3rd Qu.:10.4716   3rd Qu.:43.382   3rd Qu.:10.30   3rd Qu.:  1.60176  
 Max.   :15.6127   Max.   :80.832   Max.   :36.00   Max.   : 29.36110  
     zscore                fold      
 Min.   :-2.0312467   Min.   : 1.00  
 1st Qu.:-0.5694160   1st Qu.: 8.25  
 Median :-0.1237453   Median :15.50  
 Mean   : 0.0004815   Mean   :15.50  
 3rd Qu.: 0.3209767   3rd Qu.:22.75  
 Max.   : 4.1304203   Max.   :30.00  

plot(cv.lin$var1.pred ~ cv.lin$observed, cex = 1.2, lwd = 2)
abline(0, 1, col = "lightgrey", lwd = 2)
lm_lin <- lm(cv.lin$var1.pred ~ cv.lin$observed)
abline(lm_lin, col = "red", lwd = 2)

r2_lin = summary(lm_lin)$r.squared
rmse_lin = hydroGOF::rmse(cv.lin$var1.pred, cv.lin$observed)

## 2) Modelo exponencial
cv.exp <- krige.cv(acumulado_24h ~ x + y, locations = pluvio, model = exp.fit)
summary(cv.exp)

Object of class SpatialPointsDataFrame
Coordinates:
        min     max
x  634918.7  687966
y 7450222.8 7475508
Is projected: NA 
proj4string : [NA]
Number of points: 30
Data attributes:
   var1.pred           var1.var         observed        residual       
 Min.   : 0.06051   Min.   : 12.38   Min.   : 0.00   Min.   :-15.0221  
 1st Qu.: 3.53749   1st Qu.: 18.23   1st Qu.: 2.20   1st Qu.: -2.5773  
 Median : 5.69650   Median : 24.71   Median : 5.00   Median : -0.5739  
 Mean   : 7.01900   Mean   : 31.82   Mean   : 7.34   Mean   :  0.3210  
 3rd Qu.: 9.60323   3rd Qu.: 38.94   3rd Qu.:10.30   3rd Qu.:  2.6853  
 Max.   :20.22214   Max.   :107.38   Max.   :36.00   Max.   : 30.4535  
     zscore              fold      
 Min.   :-2.50740   Min.   : 1.00  
 1st Qu.:-0.46053   1st Qu.: 8.25  
 Median :-0.11106   Median :15.50  
 Mean   : 0.02631   Mean   :15.50  
 3rd Qu.: 0.45814   3rd Qu.:22.75  
 Max.   : 4.65631   Max.   :30.00  

plot(cv.exp$var1.pred ~ cv.exp$observed, cex = 1.2, lwd = 2)
abline(0, 1, col = "lightgrey", lwd = 2)
lm_exp <- lm(cv.exp$var1.pred ~ cv.exp$observed)
abline(lm_exp, col = "red", lwd = 2)

r2_exp = summary(lm_exp)$r.squared
rmse_exp = hydroGOF::rmse(cv.exp$var1.pred, cv.exp$observed)

## 3) Modelo Gaussiano
cv.gau <- krige.cv(acumulado_24h ~ x + y, locations = pluvio, model = gau.fit)
summary(cv.gau)

Object of class SpatialPointsDataFrame
Coordinates:
        min     max
x  634918.7  687966
y 7450222.8 7475508
Is projected: NA 
proj4string : [NA]
Number of points: 30
Data attributes:
   var1.pred         var1.var         observed        residual        
 Min.   : 0.471   Min.   : 12.83   Min.   : 0.00   Min.   :-15.95270  
 1st Qu.: 2.869   1st Qu.: 14.63   1st Qu.: 2.20   1st Qu.: -2.33169  
 Median : 6.791   Median : 19.70   Median : 5.00   Median : -0.46101  
 Mean   : 7.297   Mean   : 29.96   Mean   : 7.34   Mean   :  0.04298  
 3rd Qu.: 9.971   3rd Qu.: 40.45   3rd Qu.:10.30   3rd Qu.:  1.88845  
 Max.   :21.153   Max.   :105.30   Max.   :36.00   Max.   : 30.93771  
     zscore               fold      
 Min.   :-2.765848   Min.   : 1.00  
 1st Qu.:-0.464753   1st Qu.: 8.25  
 Median :-0.101036   Median :15.50  
 Mean   : 0.004199   Mean   :15.50  
 3rd Qu.: 0.248281   3rd Qu.:22.75  
 Max.   : 4.660612   Max.   :30.00  

plot(cv.gau$var1.pred ~ cv.gau$observed, cex = 1.2, lwd = 2)
abline(0, 1, col = "lightgrey", lwd = 2)
lm_gau <- lm(cv.gau$var1.pred ~ cv.gau$observed)
abline(lm_gau, col = "red", lwd = 2)

r2_gau = summary(lm_gau)$r.squared
rmse_gau = hydroGOF::rmse(cv.gau$var1.pred, cv.gau$observed)

## 4) Modelo Wave
cv.wav <- krige.cv(acumulado_24h ~ x + y, locations = pluvio, model = wav.fit)
summary(cv.wav)

Object of class SpatialPointsDataFrame
Coordinates:
        min     max
x  634918.7  687966
y 7450222.8 7475508
Is projected: NA 
proj4string : [NA]
Number of points: 30
Data attributes:
   var1.pred         var1.var        observed        residual       
 Min.   : 1.979   Min.   :14.15   Min.   : 0.00   Min.   :-11.4223  
 1st Qu.: 4.483   1st Qu.:20.05   1st Qu.: 2.20   1st Qu.: -4.7160  
 Median : 6.881   Median :27.94   Median : 5.00   Median : -0.5548  
 Mean   : 7.024   Mean   :28.33   Mean   : 7.34   Mean   :  0.3161  
 3rd Qu.: 9.391   3rd Qu.:35.76   3rd Qu.:10.30   3rd Qu.:  2.8856  
 Max.   :13.217   Max.   :48.97   Max.   :36.00   Max.   : 31.7935  
     zscore              fold      
 Min.   :-1.94021   Min.   : 1.00  
 1st Qu.:-0.93925   1st Qu.: 8.25  
 Median :-0.13349   Median :15.50  
 Mean   : 0.03151   Mean   :15.50  
 3rd Qu.: 0.52837   3rd Qu.:22.75  
 Max.   : 5.09689   Max.   :30.00  

plot(cv.wav$var1.pred ~ cv.wav$observed, cex = 1.2, lwd = 2)
abline(0, 1, col = "lightgrey", lwd = 2)
lm_wav <- lm(cv.wav$var1.pred ~ cv.wav$observed)
abline(lm_wav, col = "red", lwd = 2)

r2_wav = summary(lm_wav)$r.squared
rmse_wav = hydroGOF::rmse(cv.wav$var1.pred, cv.wav$observed)

# Criando uma tabela das estatística de R2 e RMSE
df.r2 <- data.frame(r2_lin, r2_exp, r2_gau, r2_wav)
df.rmse <- data.frame(rmse_lin, rmse_exp, rmse_gau, rmse_wav)
tabela <- data.frame(cbind(t(df.r2), t(df.rmse)))
colnames(tabela) <- c("R2", "RMSE")
rnames <- gsub("r2_", "", rownames(tabela))  # remove o prefixo r2 dos nomes das linhas
rownames(tabela) <- rnames  # substitui o nome das linhas simplificadas na tab original
tabela

              R2     RMSE
lin 0.1482468633 6.882779
exp 0.0765441167 7.496573
gau 0.0680566292 7.723425
wav 0.0001772582 7.991183

Criando os grids do contorno da cidade do Rio de Janeiro para intermpolação

# Importando o contorno do Rio
contorno.rio <- shapefile("dados/chuva_rio/MUNIC_2K_2022_IPP_SIRGAS.shp")

# Criando grade para interpolacao com a resolucao de 50m
r <- raster(contorno.rio, res = 50)

# Criando um objeto formato raster
rp <- rasterize(contorno.rio, r, 0)

# Trasnsformando o objeto raster no formato SpatialPixelsDataFrame
grid <- as(rp, "SpatialPixelsDataFrame")
sp::plot(grid)

Krigagem

# Colocando os dados de chuva e o grid na mesma projecao
sp::proj4string(pluvio) = CRS(proj4string(contorno.rio))

mapa_chuva_lin <- krige(acumulado_24h ~ 1, pluvio, grid, model = lin.fit)

[using ordinary kriging]

mapa_chuva_exp <- krige(acumulado_24h ~ 1, pluvio, grid, model = exp.fit)

[using ordinary kriging]

mapa_chuva_gau <- krige(acumulado_24h ~ 1, pluvio, grid, model = gau.fit)

[using ordinary kriging]

mapa_chuva_wav <- krige(acumulado_24h ~ 1, pluvio, grid, model = wav.fit)

[using ordinary kriging]

plot.lin <- plot(mapa_chuva_lin, main = "Modelo Linear")

plot.exp <- plot(mapa_chuva_exp, main = "Modelo Exponencial")

plot.gau <- plot(mapa_chuva_gau, main = "Modelo Gaussiano")

plot.wav <- plot(mapa_chuva_wav, main = "Modelo Wave")

Auto Krige

# Modelando
auto.krige = autoKrige(acumulado_24h ~ x + y, pluvio, grid, model = "Exp")

[using universal kriging]

summary(auto.krige)

krige_output:
Object of class SpatialPixelsDataFrame
Coordinates:
        min       max
x  623533.4  695333.4
y 7446574.5 7483374.5
Is projected: TRUE 
proj4string :
[+proj=utm +zone=23 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
+no_defs]
Number of points: 481750
Grid attributes:
   cellcentre.offset cellsize cells.dim
s1          623558.4       50      1436
s2         7446599.5       50       736
Data attributes:
   var1.pred           var1.var           var1.stdev     
 Min.   :-0.01083   Min.   :  0.07308   Min.   : 0.2703  
 1st Qu.: 1.50803   1st Qu.: 14.16177   1st Qu.: 3.7632  
 Median : 5.25560   Median : 22.76226   Median : 4.7710  
 Mean   : 6.66402   Mean   : 24.90904   Mean   : 4.7804  
 3rd Qu.:11.20533   3rd Qu.: 31.22818   3rd Qu.: 5.5882  
 Max.   :35.96625   Max.   :148.38326   Max.   :12.1813  

exp_var:
  np      dist    gamma dir.hor dir.ver   id
1  8  2254.389 11.19978       0       0 var1
2 27  4384.026 13.98185       0       0 var1
3 35  6073.421 36.23544       0       0 var1
4 54  8869.824 53.69680       0       0 var1
5 50 12016.434 41.83955       0       0 var1
6 47 15019.442 52.23264       0       0 var1
7 52 18659.164 85.79906       0       0 var1

var_model:
  model    psill    range
1   Nug   0.0000     0.00
2   Exp 176.6964 33077.08
Sums of squares betw. var. model and sample var.[1] 0.0003263469

# Validação cruzada
auto.krige.cv <- autoKrige.cv(acumulado_24h ~ x + y, pluvio, model = "Exp")
summary(auto.krige.cv)

            [,1]   
mean_error  0.36   
me_mean     0.04904
MAE         4.544  
MSE         57.79  
MSNE        1.789  
cor_obspred 0.2754 
cor_predres -0.3772
RMSE        7.602  
RMSE_sd     1.022  
URMSE       7.594  
iqr         4.636  

Convertendo para o formato raster - auto krige

raster_chuva <- raster(auto.krige$krige_output)
plot(raster_chuva)

• Caso queira salvar a imagem raster em um arquivo formato geotiff para ler em algum SIG por exemplo:

# Exportando o objeto da imagem
writeRaster(raster_chuva, filename = "dados/chuva_rio/chuva_auto.tiff", format = "GTiff",
    overwrite = T)

# Importando de volta para o R
raster_chuva <- raster("dados/chuva_rio/chuva_auto.tiff")

Fazendo o mapa interativo com as estações

# Importando os dados referente as estatções pluviométricas
estacoes.sf <- read_sf("dados/chuva_rio/Estac_C3_B5es_Alerta_Rio.shp")

# Convertendo UTM para Lat Long das estacoes
estacoes.longlat <- st_transform(estacoes.sf, "+proj=longlat +ellps=WGS84 +datum=WGS84")
estacoes.longlat$coords <- st_coordinates(estacoes.longlat)
estacoes.longlat$X <- estacoes.longlat$coords[, 1]
estacoes.longlat$Y <- estacoes.longlat$coords[, 2]

# Importando a malha de bairros
bairros.sf <- read_sf("dados/chuva_rio/BAIRROS_2K_2022_IPP_SIRGAS.shp")

# Convertendo UTM para Lat Long a malha dos bairros
bairros.longlat <- st_transform(bairros.sf, "+proj=longlat +ellps=WGS84 +datum=WGS84")

# Convertendo o raster da chuva para lat long
raster_chuva_longlat <- projectRaster(raster_chuva, crs = CRS("+proj=longlat +datum=WGS84"))

# Definindo Paleta de cores da superfície interpolada da chuva
pal <- colorNumeric(c("#000066", "#00c8f8", "#F0E68C", "#FFFF00", "#FF8C00"), values(raster_chuva_longlat),
    na.color = "transparent", reverse = T)

# Construindo o mapa interativo via leaflet

leaflet(data = estacoes.longlat, options = leafletOptions(attributionControl = FALSE)) |>
    # addTiles() |>
addProviderTiles("CartoDB.Positron", group = "Ruas") |>
    addProviderTiles("Esri.WorldImagery", options = providerTileOptions(opacity = 0.7),
        group = "Satélite") |>
    addProviderTiles(providers$CartoDB.Voyager, group = "Voyager") |>
    addProviderTiles(providers$Stamen.Toner, group = "Toner") |>
    setView(lng = -43.42, lat = -22.9, zoom = 10.4) |>
    # addProviderTiles(providers$CartoDB.Voyager) |>
addMarkers(~X, ~Y, popup = ~as.character(est), label = ~as.character(est), group = "Estações") |>
    ############## Polígonos dos Bairros ################
addPolygons(data = bairros.longlat, weight = 3, color = "darkblue", smoothFactor = 1,
    fill = FALSE, labelOptions = labelOptions(style = list(`font-weight` = "normal",
        padding = "3px 8px"), textsize = "13px", direction = "auto"), group = "Bairros") |>
    ########## Adicionando o raster #########################
addRasterImage(raster_chuva_longlat, colors = pal, opacity = 0.8, group = "Chuva: 1 semana") %>%
    leaflet::addLegend(pal = pal, values = values(raster_chuva_longlat), title = "Chuva Acumulada - 1 semana",
        group = "Chuva: 1 semana") |>
    ############## Controle das layers (botoes) ################
addLayersControl(baseGroups = c("Voyager", "Ruas", "Satélite", "Toner"), overlayGroups = c("Estações",
    "Bairros", "Chuva: 1 semana"), options = layersControlOptions(collapsed = FALSE),
    position = "bottomleft") |>
    ########## Desabilitando os grupos ################
hideGroup(group = c("Bairros"))